Comparison of two commercial methods for smooth-shelled mussels (Mytilus spp.) species identification

Highlights • A comparison of the sequence analysis of the H1C gene and HRM analysis of the PAP gene for species identification in the Mytilus genus was performed.• Sequence analysis of the H1C gene affects the performance of species identification.• Sequence analysis H1C reveals hybrids not detected by the HRM-PAPM.• Both methods show perfect agreement when putative hybrids are excluded from the analysis.


Introduction
Smooth-shelled mussels (Mytilus spp.) includes, at least, seven species with statistics in FAO 2019 database. Among these edible mussels, 41.5 %, 35.7 % and 14.6 % of the production came from Mytilus chilensis, M. edulis and M. galloprovincialis, respectively (FAO, 2021). It is well known that where more than one Mytilus species coexists in the same geographic area, they hybridize posing challenges to species identification (SI) (Coustau, Renaud, & Delay, 1991;Michalek, Ventura, & Sanders, 2016). Seafood international trade increased requirements regarding quality, safety, authenticity and sustainability production in food standards and regulations to give confidence to regulators, food traders and consumers, and to enable them to make informed choices (Codex Alimentarius, 2020;EU No. 1379EU No. 1380EU No. 1224/2009). These regulations require traceability "from ocean to fork", that means among others, to identify the species to which the individuals belong (Ogden et al, 2008). For example, regulation EU No.1379/2013 establishes labelling requirements for marketed seafood and foodstuffs, such as declaring the commercial designation and scientific name of the species in the product labels (D'Amico, Armani, Gianfaldoni, & Guidi, 2016;Tinacci, Giusti, Guardone, Luisi, & Armani, 2019). To enforce this regulation, methods based on DNA analysis have been employed to investigate intentional and unintentional species substitution in seafood (Fernández-Tajes et al., 2011;Santaclara et al., 2006;Verrez-Bagnis et al., 2018).
The most popular approach to specimen identification is sequencing analysis, which can be performed by different methods. One of them is the Forensically Informative Nucleotide Sequencing (FINS), which involves the estimation of sequence similarity among specimens by phylogenetic methods (i.e. Kimura-2, Tamura-Nei or Jukes-Cantor) based on genetic distances and drawing a phylogenetic tree using UPGMA or Neighbor-Joining (NJ) algorithms (Bartlett & Davidson, 1992). However, methodologies based on tree topologies perform poorly for specimen identification, therefore alternative approaches based on direct sequence comparison and genetic distances have been proposed (Collins & Cruickshank, 2013). One of these approaches is the automatic barcoding gap discovery (ABGD) for primary species delimitation (Puillandre, Lambert, Brouillet, & Achaz, 2012). ABGD begins detecting the first significant gap and uses it to partition the data. Next, limit and gap detection inferences are recursively applied to previously obtained groups to get finer partitions until no further partitioning is possible (Puillandre et al., 2012). Another approach is the "best close match" (BCM), where a specimen is successfully identified if its sequence shows the smallest genetic distance to all conspecific sequences, and it is within the 95 th percentile of all intraspecific distances (Meier, Shiyang, Vaidya, & Ng, 2006). This strategy requires determining the threshold similarity value that defines how similar a barcode match needs to be before it can be identified (Meier et al., 2006). Finally, the "all barcodes" (AB) method is a more rigorous application of the best close match strategy (Meier et al., 2006). It uses information from all conspecific barcodes in the database, instead of just focusing on the most similar ones. The barcodes are sorted by similarity to the unknown sample using the same threshold as for the best close match. The identification is achieved if at least two sequences of the query species are available in the database and when all conspecific barcodes topped the list of the best matches. This method is more confident about assigning a species name to a query in cases where multiple species names are found on the list of best matches (Meier et al., 2006). Specifically, for the genetic identification of the three Mytilus species found in Europe (M. edulis, M. trossulus and M. galloprovincialis), Inoue, Waite, Matsuoka, Odo, and Harayama (1995) developed a PCR-length polymorphism analysis from the nonrepetitive region in the polyphenolic adhesive protein gene (PAP). This region is an attractive target as it contains SNP and size polymorphisms between species within the Mytilus genus. Later, to differentiate between M. chilensis and M. galloprovincialis, Santaclara et al. (2006) added a Restriction Fragment Length Polymorphism analysis (PCR-RFLP) of the same amplicon. Always targeting this gene, Jilberto, Araneda, and Larraín (2017) developed a High-Resolution Melting (HRM) analysis (HRM-PAPM) to identify M. chilensis, M. galloprovincialis, M. edulis and their F 1 hybrids. By the HRM analysis is possible to obtain the genotypes of single nucleotide polymorphisms (SNPs) and detects length polymorphisms in small amplicons, showing distinguishable melting curves, permitting discrimination among species (Verrez-Bagnis et al., 2018). The HRM-PAPM analysis allows for distinguish, M. galloprovincialis (genotype GG, amplicon size 116 bp), M. edulis (genotype GG, 170 bp), and M. chilensis (genotype TT, 116 bp) (Jilberto et al., 2017). HRM is ideal for laboratory analysis because is fast, accurate and less expensive in comparison with other DNA based methodologies (Jilberto et al., 2017). Another method used in the identification of Mytilus species is the comparison of the partial sequence of the histone H1C gene (Eirín-López, González-Tizón, Martinez, & Méndez, 2002). Specifically, H1C gene is used to separate M. chilensis from other Mytilus species (M. galloprovincialis, M. edulis, and M. trossulus). Both, the HRM-PAPM method and the H1C sequencing are easy to apply in routine testing to identify Mytilus species in traded seafood. However, basing SI on the analysis of a single gene, known as the mono-locus approach, can give contradictory results due to differences in the evolutionary rate of the analyzed genes (Larraín, González, Pérez, & Araneda, 2019;Väinölä & Strelkov, 2011). Also, not all markers can differentiate among all species (Larraín et al., 2019). To overcome these problems related to the analysis of a single marker, the multi-locus analysis integrating genomic information from SNP panels to identify M. trossulus, M. galloprovincialis, M. edulis, and M. chilensis has been used by Larraín, Zbawicka, Araneda, Gardner, and Wenne (2018) and Wenne et al. (2022) with 49 and 54 SNPs, respectively.
In this work, we tested and compared the performance of the HRM-PAPM method and the H1C sequence analysis for species identification in smooth-shelled mussels Mytilus galloprovincialis and M. chilensis using individuals previously identified with a panel of 49 SNPs.

Mussel samples, DNA extraction and species identification
Fifty-two samples of the Mytilus mussels had obtained from six growing centers in the Reloncaví Sound area (− 41.700, − 72.833), where the Chilean mussel aquaculture industry is located and one in Peel Island (− 50.842, − 74.011). Nine fresh individuals were obtained from a growing center in the Dichato Bay (− 36.538, − 72.957), 570 km away from the aquaculture area. Reloncaví Sound and Peel Island samples were collected in 2009 and 2013 (n = 26 per year), and Dichato samples were collected in 2009 (Table S1). All samples were processed fresh up to 24 h after collection, dissected and approximately 200 mg of mantle edge tissue was fixed with 95 % ethanol and stored at − 20 • C until DNA extraction. DNA was obtained by the phenol-chloroform method, adapted for mussels (Larraín, Díaz, Lamas, Vargas, & Araneda, 2012). The DNA concentration was estimated with a NanoDrop ND-2000 spectrophotometer (Thermo Fischer Scientific). DNA integrity was assessed using agarose gel electrophoresis (0.7 % w/v). Mussel species were previously determined with a 49 SNPs panel by Larraín et al. (2018), resulting in nine as M. galloprovincialis (Mg) and 52 as M. chilensis (Mch).

High resolution melting (HRM-PAPM) analysis
The HRM-PAPM analysis was performed in an Eco Real-Time PCR System 4.0 (Illumina®) and Mic qPCR Cycler (Bio Molecular Systems). In all HRM analyses, M. galloprovincialis, M. edulis, and M. chilensis reference samples were included as controls. Also, negative control without DNA was included in all runs. Species identification was performed through clearly distinguishable melting curves, for extended protocol see Jilberto et al. (2017). The validation of the HRM-PAPM analysis was published in Quintrel et al. (2021), and a validation summary is included in Table S2.

Histone H1C gene sequencing
A partial sequence of approximately 400 bp of the H1C gene was amplified by PCR using the primer set H1CF (5 ′ -CATCATGGCCAACTT-CAACG-3 ′ ) and H1CR (5 ′ -GGCTGAATAGCCTCTGCAGA-3 ′ ) (Pérez-García, Morán, & Pasantes, 2014). The final product length was checked by electrophoresis on a 2 % agarose gel. PCR and sequencing were performed in the facilities of Eurofins Genomics GmbH (Ebersberg, Germany). All reactions were carried out in a 25-µL volume containing 1.5 U AmpliTaq Gold DNA Polymerase (Applied Biosystems®), 2.5 µL 10 × Gold Star Buffer (Promega Corporation®), 1 to 5 ng DNA and 0,2 µM of each primer). The PCR reaction was performed on a GeneAmp 9700 thermocycler (Applied Biosystems®). Initial denaturation was performed at 95 • C for 12 min; followed by 15 cycles of denaturation at 96 • C for 25 s, annealing at 55 • C for 20 s, and extension at 72 • C for 30 s; followed by 25 cycles of denaturation at 96 • C for 20 s, annealing at 58 • C for 20 s, and extension at 72 • C for 30 s with a final elongation step at 72 • C for 3 min. Double-stranded PCR amplicons sequencing was performed in an ABI 3130xl genetic analyzer (Applied Biosystems®).

Direct sequence comparison (DSC)
The raw chromatograms of the 61 individuals were manually edited in Geneious 11.5 (Biomatters) and aligned as described in FINS analysis (2.3.1). According to this alignment analysis, the species of M. chilensis and M. galloprovincialis individuals were assigned considering the genotype at polymorphic sites (single or double pick in chromatograms) shared between both species. In this case, individuals heterozygous showing double picks in these sites were considered as putative hybrids. Therefore, two datasets were defined: the first one (dataset1) containing sequences from all the 61 individuals of the Mytilus genus sampled in this study and described in 2.1. The second dataset (dataset2) contains the abovementioned sequences, excluding the 12 putatively hybrid individuals determined from the analysis of the raw chromatograms (see Table 2 and results in 3.2.2). An individual was assigned to a species when the percentage of similarity of its H1C sequence was greater than 98 % compared with the other conspecific sequences (Armani et al., 2015;Barbuto et al., 2010;Hebert, Ratnasingham, & DeWaard, 2003). A summary of the validation for the DSC analysis on the H1C gene is included in Table S2.

Other sequence analysis methods
Also, we tested the other three methods commonly used in sequence comparison analysis to identify M. chilensis and M. galloprovincialis individuals. First, the ABGD method was run 4 times from its webpage: htt ps://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html with both evolution models (Jukes-Cantor JC69 and Kimura 2-parameter K80) and a Prior Intraspecific divergence value between 0.1 and 0.001 (Puillandre et al., 2012). The second and third methods were the BCM and the AB, both applied to estimate a threshold similarity value to consider individuals as conspecific (Meier et al., 2006)¡. They were run twice in the TaxonDNA/Species identifier ver. 1.8 software with Kimura 2-parameters correction using the two previously described datasets (Meier et al., 2006).

Agreement between the HRM-PAPM and histone H1C sequence analysis
The classification agreement between the HRM-PAPM and each of the Histone H1C sequence analyses described in 2.3 (DSC, BCM, AB and ABGD methods), was evaluated using individuals from dataset1 and dataset2. We calculate Cohen's Kappa coefficient (κ) and the Matthews correlation coefficient (MCC) to evaluate the quality of SI obtained by the different molecular assays. κ statistic measures the agreement between two methods that classify items into mutually exclusive categories (Rotondi, 2018), and it was estimated with the 95 % confidence interval (Cohen, 1960), using the "fmsb" R package version 0.6.3 (Nakasawa, 2018). FDR correction was applied to κ p-values to avoid the alphainflation produced by multiple testing (García, 2004). MCC is a measure of the quality of classification agreement between two categorical variables (Matthews, 1975) that was estimated using the "mltools" R package version 0.3.5.

HRM-PAPM analysis
This method based on the melting curves from the 61 samples, classified nine individuals as M. galloprovincialis and the remaining 52 as M. chilensis, without any evidence of hybridization ( Figure S1).

Fins
The phylogenetic tree shows two principal clades (Fig. 1), the first one containing only two individuals from data uploaded by Eirín-López et al. (2002) on GenBank as M californianus (AJ416421), and M. chilensis (AJ416422). The second clade grouped the remaining 75 individuals, among them, only M. trossulus was separated significantly with H1C gene analysis. The separations within this clade did not reach the minimum threshold confidence to be considered significant.

DSC
After aligning the H1C 397 bp amplicon from all sequences available from GenBank and this work, 49 polymorphic sites were found. Despite    these polymorphic sites, all sequences showed at least 91.44 % of similarity (Tables S3 and S4). The lowest similarity value between M. galloprovincialis and either M. edulis (download from GenBank) or M. chilensis was 97.73 %. (Table S4). Besides, in some cases, sequences corresponding to M. edulis and M. galloprovincialis were identical (Table S3).
Twelve individuals (GenBank accession numbers MT949820, MT949828 to MT949837 and MT949795) showed chromatograms curves with double peaks at five sites. These sites were: 115,151,193,259,322 (Table 2). In M. galloprovincialis these sites corresponded to T, A, C, G and A respectively, whereas in M. chilensis were C, G, T, A and C, respectively. The 12 individuals that presented double peaks in these five sites, showed the two nucleotides corresponding to the polymorphisms found between M. chilensis and M. galloprovincialis mentioned above, and were considered heterozygous for these sites and putative hybrids.
The sequence matrix showed a 97.73 % to 98.74 % of similarity between M. chilensis and M. galloprovincialis (Table S3 and S4). Considering the twelve sequences of putative hybrids carrying double picks in the five sites, the similarity with M. chilensis sequences ranged from 98.36 % to 99.87 % and with M. galloprovincialis ranged from 97.86 % to 99.37 % (Tables S3 and S4). Therefore, pairwise similarities among individuals from both species and their putative hybrids were higher than 97.73 % (Table S3 and S4).
Using dataset2 (excluding individuals whose sequences showed double picks) DSC analysis successfully assigned all samples to the correct species when a sequence similarity threshold over 99 % was considered.

Other sequence analysis methods
The ABGD method was not able to separate the individuals from dataset1 into groups according to the nominal species classification, the prior intraspecific distance ranged between 0.0028 and 0.0077. The BCM and the AB methods calculated a threshold similarity value of 1.28 %. BCM analysis correctly recognized the species from 55 sequences (90.16 %), four (6.56 %) were ambiguous (MT949829, MT949830, MT949832 and MT949837) and two (3.28 %) were incorrectly classified (MT949828 and MT949836). The ambiguous or incorrectly classified sequences corresponded to putative hybrids found in the DSC analysis. However, six sequences of those individuals who presented double picks (MT949795, MT949820, MT949831, MT949833, MT949834 and MT949835) were classified as M. chilensis. The AB method classified all the 61 individuals from dataset1 as ambiguous (Table 3).
Using dataset2, the ABGD method detected a gap in the intraspecific and interspecific distances in H1C sequences between M. chilensis and M. galloprovincialis. The BCM and the AB methods calculated a threshold similarity value of 0.25 %, to identify barcodes as conspecific. The BCM analysis successfully assigned to the species, all the 49 individuals. However, three sequences matched with the species but with interspecific distances outside the 0.25 % threshold (MT949794, MT949825 and MT949827). Besides, the AB method was able to correctly identify the seven M. galloprovincialis and 39 M. chilensis, while three individuals were classified as "no-match". These correspond to the same three individuals successfully identified by the BCM method outside the threshold.

Discussion
Molecular taxonomic identification or specimen assignment, i.e., to assign an individual to a species, is important for food safety and authenticity, wildlife forensic, conservation, property and consumers rights protection. In the food sector including, regulators, processors, retailers and consumers, there is special concern about species identification (Armani et al., 2015;Verrez-Bagnis et al., 2018). It is widely accepted the need for proper methods to support law enforcement, that could be easily implemented by private or government laboratories. However, when different analytical methods are used for the same purpose (in this case, SI) the results may not be concordant among them. Therefore, harmonization is necessary to avoid stakeholders being affected by the discrepancies among the different analytical methods (Coleman & Fontana, 2010). The harmonization among commercial methods currently used for mussel SI begins with a comparison of their performance. Discrepancies among methods could arise from the different molecular markers analyzed, but also from the sequence analysis methodology employed. In this study, we compare two monolocus commercial methods used to identify the most traded Mytilus mussel species: the HRM-PAPM and the sequence analysis of the H1C barcode gene.
The high genetic similarity among Mytilus spp. has been revealed by many markers used for DNA barcoding (COI, 5S rDNA, ITS-1, ITS-2 and NTS) (Giusti et al., 2020;Santaclara et al., 2006;Tinacci et al., 2018), which were not able to separate the different species. However, the polyphenolic adhesive protein (PAPM) and the histone H1C genes, have been shown to successfully differentiate Mytilus mussels (e.g., Santaclara et al., 2006, Fernández-Tajes et al., 2011, Eirín-López et al., 2002. SI results obtained by the HRM-PAPM and H1C gene sequence analysis showed a slight to a substantial agreement among them, using the complete dataset (dataset 1). Discrepancies between mono-locus methods are expected because they target genomic regions with different evolutionary histories (Larraín et al., 2012;Väinölä & Strelkov, 2011). The PAPM is a single copy gene, while the H1 genes (including H1C) have evolved by gene duplication with close to one hundred copies per haploid genome in M. galloprovincialis (Eirín-López et al., 2002) and M. edulis (Drabent et al., 1999). Besides, in contrast to the highly conserved family of core histones (H2A, H2B, H3, and H4), the H1 histones are less conserved during evolution (Drabent et al., 1999).
H1C sequences showed very high similitude among the Mytilus taxa analyzed, therefore, there was no barcoding gap in the H1C gene to separate intra-from inter-specific sequence variation when considering all the individuals analyzed. Probably, this is the reason why FINS methodology based on the H1C gene was not useful for SI in smooth- shelled mussels (Fig. 1). Also, considering all the problems of tree-based identification techniques discussed by Meier et al. (2006) and Collins and Cruickshank (2013), FINS is not a recommendable approach to perform traceability in Mytilus species. On the other hand, the DSC method uses a similarity criterion relying on an arbitrary threshold to consider two sequences belonging to the same species. Usually, a 1 or 2 % divergence is a reasonable rule-ofthumb in most cases, but it is prone to produce different rates of falsepositive and false-negative (Collins & Cruickshank, 2013). In this study, a 1 % divergence threshold in pairwise similarities between M. chilensis and M. galloprovincialis H1C sequences, was suitable to correctly identify all samples, excluding putative hybrids (dataset2). The ABGD method was not able to identify a barcoding gap between M. chilensis and M. galloprovincialis sequences when all individuals were considered, as is in a real scenario in mussel traceability. This sequence analysis method only worked when putatively hybrids were previously excluded from the analysis. As expected, the BCM and AB methods showed a broader threshold similarity value (1.28 %) when all individuals were considered, compared with the 0.25 % threshold obtained when putatively hybrids were excluded. When these two methods were applied to identify the species in all individuals, the performance was lower. Moreover, the strictest criteria of the AB method, classified all individuals as ambiguous. It is important to consider that using the same dataset, sequence analysis methods could give different results. Therefore, the election of the sequence analysis method for SI must be considered during the standardization and validation process.
The sequences for the H1C gen in Mytilus spp. available in the Gen-Bank, were published in two papers. Eirín-López et al. (2002) uploaded sequences from M. edulis, M. galloprovincialis, M. trossulus, M. chilensis and M. californianus, whereas Drabent et al. (1999) contributed with nine M. edulis sequences for this gene. Until now, the NCBI Reference Sequence Database contains no sequences for the H1C gen from Mytilus spp. The reference sequences have their provenance and validity reviewed and checked and their GenBank ID begins with NC. This quality check is extremely important to avoid taxonomic misidentifications. Reviewing the H1C sequences published in GenBank, we realized that the sequence AJ416422 uploaded as M. chilensis by Eirín-López et al. (2002) grouped with the M. californianus sequence and not with M. chilensis (Fig. 1). Besides, 27 polymorphic sites were present between our sequences of M. chilensis and the sequence AJ416422, indicating a possible error in the species assigned to the last sequence in GenBank. Mistakes uploading sequences to GenBank, especially at the species level (but not at the genus level) are common, affecting the taxonomic reliability of this database (Leray, Knowlton, Ho, Nguyen, & Machida, 2019). Both species are very different in shell morphology, M. californianus is a ribbed mussel, while M. chilensis is a smooth shelled mussel. Unfortunately, the specimen used by Erin-Lopez do not have morphological data linked to the GenBank record. The availability of Table 4 Kappa concordance, confidence interval of 95%, agreement and Matthews correlation coefficient (MCC) between Species Identification methods based on the HRM-PAPM / Me15-16 and H1C barcode gene Direct Sequence Comparison (DSC), Best close Match (BCM) all barcode (AB) and Automatic Barcode Gap Discovery (ABGD) methods for a) dataset1(all sequences) and b) dataset2 (without putative hybrids).  reference sequences obtained from vouchered specimens in cured databases is necessary to standardize sequence-based methods as traceability tools and to avoid taxonomical uncertainties. The species from all the individuals whose H1C sequences were obtained in this work and uploaded to GenBank were previously checked by the 49 SNP panel (Larraín et al., 2018) and corresponded to M. chilensis and M. galloprovincialis. Food traceability and authenticity in Mytilus mussels is challenging due to the natural hybridization observed in many areas where these species are cultured (Michalek et al., 2016). In the case of M. chilensis and M. galloprovincialis the presence of hybrids has been highlighted in molecular identification studies conducted on cultured specimens in Arauco and Reloncaví gulfs (Larraín et al. 2012(Larraín et al. , 2019) and on exported Chilean mussel products collected from retail (Giusti et al., 2022). To overcome this challenge, new molecular techniques are necessary to separate hybrids from the parental taxa in an efficient a cost-effective way. However, according to our results, mono-locus assays are not efficient for this purpose, because they gave contradictory conclusions. Methods based on a multi-locus approach (Larraín et al., 2019) or genomic technology must be developed and validated to be used by regulatory agencies and the food industry. The higher the complexity of a method, the more time it consumes, affecting the practicability required in routine analysis. Consequently, using a small panel of informative SNPs is a good option to address the presence of hybrids (Quintrel et al., 2021). For example, all individuals used in this work were previously genotyped with a 49 SNP panel without any evidence of being hybrids. However, H1C sequence analysis showed strong evidence of hybridization in twelve individuals, reflecting the complexities of hybrids identification.
The SI discrepancies between these H1C and PAP genes were mainly because the putative hybrids detected by the H1C gene sequence analysis were not classified as such by the HRM-PAPM method, although this latter can also distinguish hybrids between these species (Jilberto et al., 2017).
The HRM-PAPM method showed total agreement with the 49 SNPs panel to identify M. chilensis from M. galloprovincialis. The H1C sequence analysis is helpful to differentiate between both species when the appropriate sequence analysis method is used.. Our results show that H1C can reveal some hybrid not detected by the HRM-PAPM, finding that is highly possible because Mytilus mussel hybridizes, their genome has a size of ~ 1.6 Gb and of course we are screening only a very low fraction of it with our SI methods.

Conclusions
A high similarity of H1C gene within Mytilus spp. is evident and it hampers the differentiation between M. edulis and M. galloprovincialis. Nevertheless, the study confirms that the application of the BMC method allows both discrimination between the two species and the detection of hybrids.
Based on the H1C gene, the sequence analysis method affects SI outcomes. The FINS, ABGD and AB methods were not useful to identify M. chilensis, M. edulis and M. galloprovincialis specimens.
The different levels of agreement between the 49 SNP panel and both SI methods (HRM-PAPM and H1C gene sequence analysis) highlight the need for the standardization of molecular tools. The presence of hybrids is a realistic scenario in Mytilus mussel aquaculture that makes SI more complex. Our results indicates that analytical tools based on a multilocus approach are needed to enface the challenge of the traceability of smooth-shelled mussels.
Each assay laboratory must decide what is the best method for mussel SI depending on the equipment available and if it can easily be combined with other methods used in the laboratory. If sanger sequencing is a routine procedure, we recommend the sequencing of the H1C gene followed by analysis with the BCM method. On the other hand, if the qPCR analysis is the routine procedure HRM-PAPM is the recommended method.
Our results highlight the need for standardized molecular tools to perform SI in smooth-shelled mussels, as well as the use of a multi-locus approach.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.